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Abstract 

We present the first 2+1 flavour lattice QCD calculations of pseudoscalar flavour-singlet propaga- 
tors using improved staggered fermions. We explore the relevant techniques and discuss prospects 
for the larger scale studies now in progress. The disconnected correlator is shown to have a highly 
non-Gaussian distribution and reliable estimates of the errors require care. 
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I. INTRODUCTION 



The high mass of the rj meson relative to the light pseudoscalar mesons is thought to 
be due to topology and the complexities of the QCD vacuum 

HQ. 

The mass of the 

rf meson is a clear experimental signal of the complexities of the QCD vacuum that are 
normally obscured by confinement, so a first principles calculation of the mass would be an 
important milestone in taming non-perturbative QCD. A robust first-principles calculation 
of the spectrum of the 77/7/ system should help shed light on the workings of the fermion 
sea, and the mechanism by which quark loops elevate the mass of the singlet meson above 



that of the octet meson 



,I2|. 

The computation of the mass of the 7/ meson has always been one of the goals of the 
program of hadron spectroscopy from lattice QCD calculations. There have been many 
lattice studies that have computed the mass spectrum of pseudoscalar flavour- singlet mesons 
with N f = 2 flavours of fermions Q, 2, B, B, Hfl 1 3, 11 1 - Recently the JLQCD/CP-PACS 
collaboration reported on a preliminary calculation of the 77 and 7/ mesons with 2 + 1 flavours 
of Wilson fermions [12j . 



Improved staggered fermions have provided some of the most impressive results in hadron 
spectroscopy for flavour non-singlet quantities [3, Q. The next stage is to use improved 
staggered fermions to compute the mass spectrum of light singlet pseudoscalar mesons. Im- 
proved staggered fermions have the advantage that they are fast to simulate and allow lattice 
QCD calculations with high statistics at relatively light quark masses, both requirements 
for accurate calculation of the masses of singlet pseudoscalar mesons. 

However, lattice QCD calculations of the masses of the 77 and 7/ mesons face several 
challenges. The lattice QCD calculations of the flavour singlet mesons involve disconnected 
correlators which are computationally more expensive to compute than connected correla- 
tors. Also disconnected diagrams are inherently noisy in lattice simulations, so high statistics 
are required. The disconnected diagram for the pseudoscalar mesons is related to the topo- 
logical charge of the gauge configuration which in some cases has been seen to have longer 
autocorrelation times than other quantities in lattice QCD simulations. To be sure statis- 
tical fluctuations and any autocorrelation times are well under control high statistics are 
required for an accurate lattice QCD calculation. 

Although lattice QCD calculations that use improved staggered fermions have been suc- 
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cessfully tested against experiment for many quantities Q], the formalism has not been 
proved to be fully correct because of the "fourth-root trick" emp loyed to achieve the required 
flavour structure. There is increasing theoretical work [jj], [if]] that suggests that there are 
no problems with the improved staggered fermion formalism in the continuum limit, how- 
ever the issue has not been completely settled. The rp]' system is a good place to check the 
validity of the rooting of the sea quark determinant, because of the important role the sea 
quark loops play in raising the mass of the rj' meson. Creutz has recently claimed that the 
mass of the 77' meson is a place where the improved staggered fermion formalism may not 
work [nj, Q- Kilcup and Venkataraman [,9j have studied the flavour singlet pseudoscalar 
meson with naive staggered fermions in unquenched QCD. 

The 77 and 7/ mesons can both be studied with lattice QCD calculations that include the 
dynamics of 2 + 1 flavours of sea quarks. In a quark model treatment [19] of the 77 and 7/ 
mesons in terms of light and strange quarks: 

i] ps 0.58(m75M + djsd) — 0.57 S75S 
7/ ps 0.40(m7 5 w + d>y 5 d) + 0.82 s^s . 

This suggests that both the 77 and 7/ mesons contain mixtures of light and strange quarks. 
The uu + dd and ss interpolating operators will both couple to the 77 and 77' mesons; in the 
singlet pseudoscalar channel the lightest state will be the 77 and the first excited state the 
77' meson. This is an additional complication of lattice QCD calculations that include 2+1 
flavours of sea quarks over those which only include 2 flavours of sea quarks. In Nf = 2 
full lattice QCD calculations the emphasis is on finding the large mass splitting between the 
singlet pseudoscalar meson and the flavour non-singlet flavour mesons. With 2+1 flavours 
of sea quarks it also interesting to see how the mass of the 77 meson is reproduced because 
this is sensitive to the mixing of (775(7 and S75S loops. The reader is referred to the review 
by Feldmann [^(J for a recent discussion of 7777' mixing. The modern way to describe 7777' 
mixing is via leptonic decay constants [20|. See 4j for an attempt to compute the relevant 
decay constants in Nf = 2 unquenched QCD. 

In this paper we report on the lattice methods required to study the singlet pseudoscalar 
mesons using improved staggered fermions. This work is the necessary starting point for a 
project that uses a large number of Nf = 2 + 1 configurations. For some cross-checks on 
the unquenched calculations we also study the singlet pseudoscalar meson in quenched QCD 



and briefly investigate the role of topology. 

The general plan of the paper is as follows: In Section [Til we describe the theory behind 
computing the correlators of singlet pseudoscalar mesons using improved staggered fermions. 
Next, in Section HTl] we describe the details of the calculation and report on the algorithmic 
work to compute the required disconnected diagrams. Section IIVI contains the results from 
the quenched and unquenched calculations and in Section [V] we discuss the statistics of 
singlet and non-singlet correlators. After that we estimate the number of configurations 
required to get a given accuracy in Section [VIJ The paper ends with our conclusions in 
Section ETU 

II. THEORETICAL BACKGROUND 

For notation, we use SP for the singlet pseudoscalar meson and NP for the non-singlet 
pseudoscalar meson. The additional complications due to staggered fermions are discussed 
later in this section. We begin with the expression of the pseudoscalar singlet propagator 
for Nf degenerate flavours of fermions: 

G S p(x',x) = (^q i (rfhBqi(x')^q j (x)'y 6 q j (x)). (1) 
i=i j=i 

This propagator gives rise to two different types of diagrams. There are Nf connected 
diagrams: 

, * , 

(^Qi^'h^qiix^^qjix^qjix)), (2) 



and Nf disconnected terms: 



The connected term in the Nf = 2 flavour symmetric case is the same as that for the 
(non-singlet) pion propagator. 
So we can write: 

G SP (x, x) = N f C(x', x) - N 2 f D(x', x), (4) 

where the extra fermion loop in the disconnected diagram gives rise to the relative minus 
sign, and 

G NP {x',x) = N f C{x\x). (5) 
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A. D/C ratio 



In full Euclidean QCD we expect that the pion propagator should decay exponentially 
at large time separations as 

G NP (t) = Ae- mNpt (6) 
and similarly for the singlet propagator: 

G SP (t) = Be- mspt . (7) 

As noted above the singlet propagator contains a disconnected and connected part, the latter 
being proportional to the pion propagator: 

G SP {t) = Ae~ mNpt - N 2 f D(t) . (8) 

Taking the ratio R(t) of the disconnected to connected parts at large time-separation At 
then suggests 

NjD(t) B , 

[ ) - N f C{t) ~ 1 A 6 ' [J) 

This derivation of ([9]) requires one to assume the same action is governing both sea quarks 

and valence quarks — in systems where the physics of the valence and sea quarks differ this 

expression may no longer apply. An important case of this is in the quenched limit where 

the D/C ratio is Q 

R(t) = A' + B't . (10) 

Another important case, relevant for staggered fermions, is when the number of flavours 
Nf of sea quarks differs from the native number of valence fermions, N v = 4. In practice, 
this is the case for most staggered simulations — while the fourth-root trick reduces four 
degenerate sea flavours to the desired Nf, there remain four native valence flavours, or tastes, 
which contribute to singlet propagators. In the connected contribution (top diagram in Fig. 
[TJ a single valence loop connecting the endpoints (meson interpolating operators) has four 
tastes of fermions circulating, introducing a factor of four. However in each of the terms of 
the disconnected contribution (bottom diagram in Fig. [TJ there are two valence loops, one 
at each endpoint, each contributing a factor of four. The D/C ratio therefore naively has 
an extra factor of N v and, if the staggered formulation (incorporating a fourth-root of the 
determinant) correctly reproduces Nf sea flavours, we expect: 

N 2 f D(t) 

Kit) = ihrri = ^,(1 - Ke-<- msp ~ m » p ») . (11) 



In the numerical work described in the following sections we implicitly rescale the discon- 
nected contributions D by a factor of 1/4 so as to correct for the extra valence tastes, a 



procedure which Sharpe calls "valence rooting" [15] 



and was used by Venkataraman and Kil- 



cup in |9|. We will continue to refer to use the NP and SP notation, rather than introduce 
some new acronyms based on taste singlet notation. 

There remains a possibility that the fourth-root trick could introduce some other patholo- 



171 ]. e.g. the wrong number of sea flavours or a mismatch between 



211 ] . Such issues should, in principle, be detectable in the 



gies at finite lattice spacing 
sea and valence fermion masses 
D/C ratio. 

For the case of Nf degenerate flavours there is only one possible interpolating operator 
for flavour singlet mesons for non-degenerate flavours of sea quarks there are more possible 



interpolating operators. In quark model inspired discussions of the 7777' mixing 19|, the 
calculations use the SU(3) octet rj 8 and singlet rj basis states 

M75M + d^d + S75S 



m 

Vs 



U'jsU + d^ 5 d — 2s7 5 s 
V6 



(12) 
(13) 



These basis states mix, because the SU(3) symmetry is broken by the large mass of the 
strange quark, to form the physical 77 and rf mesons. 

Another possible basis for flavour singlet pseudoscalar mesons is the quark basis, [21 

M7 5 m + d<y 5 d 



V2 



S75S. 



(14) 



The r] q and rj s will both also couple to the 77 and rj' mesons, and are valid interpolating 
operators for flavour singlet pseudoscalar mesons. Each of these interpolating operators (77^, 
r) s , ?7o and r/g) will have the 77 meson as the ground state and the rj' as the first excited state. 

In principle the _+ glueball could mix with the 77 and 77' mesons. This would mean 
that a 0~ + glueball interpolating operator should also be used to study 77 and 77' mesons. 
However quenched glueball simulations [22j, [23] suggest the mass of the _+ glueball to 
be around 2.6 GeV. Since this is far from the mass of the 77', we do not consider the 0~ + 



glueball interpolating operator in what follows. Hart and Teper 



24| studied the _+ glueball 



<r 3 C 3> 



FIG. 1: The valence contribution to the connected (top) and disconnected (bottom) diagrams. 
There are natively N v = 4 valence tastes in the staggered formulation, so each of these loops 
introduces a factor of four. 

correlators in unquenched QCD and, despite large statistical errors, claimed the results were 
similar to those from quenched QCD. 

The rj does not decay via the strong interaction, hence the mass of the rj is a gold-plated 
quantity [3] that should be possible to compute very accurately. The rf decays via the 



strong interaction but it has a small width of 0.2 MeV [25j (relative to that of the p meson, 
for example). The dominant strong interaction decay of the rj is via the decay to rjirir and 
thus we expect that this decay channel is not open for the masses we use. Thence we expect 
that an accurate calculation of the rj' mass should be possible in principle. 

For our Nf = 2 + 1 simulations, if we use the rjo interpolating operator, we must generalise 
flU) and (jHJ) to the case of non-degenerate flavours. The full propagator becomes 

Gspsusit) = 2C qq (t) + C ss (t) - AD qq (t) + AD qs (t) + D ss (t), (15) 

and likewise the D/C ratio is 

Rm AD qq (t)+AD qs (t)+D ss (t) 

R{t)sU3 ~ 2Cjt) + C ss (t) ' (16) 

where D qq , D qs , and D ss represent the disconnected correlators constructed with two light 
quark loops, one light and one strange quark loop, and two strange quark loops, respectively. 
Likewise C qq and C ss are the connected correlators measured with light and with strange 
quark masses, respectively. 



7 



B. Fitting the singlet propagator 



With 2 + 1 flavours of light quarks we have the potential to investigate the spectrum of 
both the 7] and the rj' states. As discussed in Section [U in 2 + 1 flavour calculations, the rf 
meson is the first excited state contributing to the pseudoscalar correlator. One approach 
would be to measure the appropriate connected and disconnected correlators, assemble the 
full singlet propagator in ffl5l) and then fit to a multi-exponential form: 



G(t) SP = A ie - mspit + A 2 e- msp2t . 



(17) 



However, fitting a multi-exponential expression with several parameters to data consisting of 
a single correlator can be difficult when that data is precise, and even more so when dealing 
with inherently noisy disconnected correlators. 

In a multi-channel approach, treating the non-strange (q) and strange (s) contributions 
separately, using the quark basis (HM one would form the matrix 



G(At) 



and then fit 



vt(At) Vq (0) vt(&t) Vs (0) 

C qq (At) - 2T> qq (At) -V2B qs {At) 
-V2V sq (At) C ss (At) - T> ss (At) 



(18) 



G(t) = A T e~ mt A (19) 

where A is a matrix of amplitudes and e~ mt is a diagonal matrix. This variational approach 
to fitting is one of the main methods of extracting masses and decay constants of excited 
mesons from lattice QCD calculations. The variational method is reviewed in 26] and used 



in 



27 



28 



29 



30 



3l\ . The recent work on the rj and rj' mesons by the CP-PACS/JLQCD 



collaboration used the variational method 



12|. 



Further improvement might be obtained by replacing the Cs and Ds in ffl8|) by matrices 
of correlators formed with the various combinations of fuzzed and point source and sink 
operators. We leave investigation of such variational fits for future work. 
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N, 


w/g 2 


L 3 x T 






am val 


^cfg 







8.0 


16 3 x 32 







0.020 


76 




2 


7.2 


16 3 x 32 




0.020 


0.020 


268 







8.00 


20 3 x 64 






0.020 


408 







8.00 


20 3 x 64 






0.050 


408 — 


-> 6154 


2 


7.20 


20 3 x 64 




0.020 


0.020 


547 




2+1 


6.76 


20 3 x 64 




0.007, 0.05 


0.007, 0.05 


422 




2+1 


6.76 


20 3 x 64 




0.010, 0.05 


0.010, 0.05 


644 




2+1 


6.85 


20 3 x 64 




0.05, 0.05 


0.05, 0.05 


369 




TABLE I: Ensembles used for singlet calculations. 


Label 


MILC 


SO SI 


S2 


S3 S4 


S5 S6 S7 


S8 S9 


all 


Configs 


408 


2455 413 


414 


409 411 


409 410 409 


412 412 


6154 



TABLE II: Quenched configurations in the MILC and extended ensemble streams. 
III. SIMULATION AND MEASUREMENT 
A. Configuration ensembles 

We performed singlet correlator measurements on the Nf = 0, 2, and 2 + 1 ensembles 
listed in Table [H We used 16 3 x 32 lattices for algorithm tuning as discussed below, and 
20 3 x 64 lattices for physics measurements. The latter are primarily the so-called "coarse" 
lattice gauge configurations generated by the MILC collaboration |32j with the "Asqtad" 



33 



34 



35 



36]. 



improved staggered action 

For analysis of the fluctuations in disconnected correlators (Section Ej) we extended the 
(3 = 8.00 quenched ensemble from 408 configurations to 6154 configurations. From the 
final configuration in the MILC ensemble, configuration 4090, we initiated ten separate 
Markov streams with 10 sweeps, each consisting of four over-relaxation steps and and a 
quasi-heatbath hit, between saved configurations. These are listed in Table ILT1 



9 



B. Connected and disconnected correlators 



The numerical calculation of the singlet propagator is performed in two parts — the 
calculation of the connected contributions, and the calculation of the disconnected contribu- 
tions. The former is relatively straightforward — inversions on point source vectors produce 
quark propagators which are then multiplied together with an appropriate meson operator 
to produce a meson propagator. 

There are in principle two different meson operators which couple to the flavour-sin glet 
pseudoscalar meson. These are the (75 ® 1) and the (7475 <8> 1), using the Kluberg-Stern [37] 
notation for the staggered meson operators. The former has the quark and anti-quark sources 
separated by four links, on opposite corners of the hypercube, while the latter has the quark 
and anti-quarks separated by three links, sited at opposite corners of the spatial cube. We 
write the operator as (75 (g> 1) to show that the meson has pseudoscalar Dirac structure, but 
is a singlet in Kogut-Susskind taste space. 

As is generally the case with staggered meson operators, the operator which couples to 
the (7475 <S>1) meson also couples to a parity partner state, which is in this case the (107475) 
scalar meson. The parity partner of the (75 ® 1) however is exotic and therefore makes no 
contribution to the (75 ® 1) propagator measured on the lattice. Furthermore, a variance 
reduction trick we discuss in Section IIII CI applies only to the (75 (g> 1) operator. For these 
reasons we use the (75 ® 1) state exclusively in this work. 

To apply the A 76< g>i operator to a source vector or fermion propagator, we covariantly 
and symmetrically shift the source by one lattice unit in each of the four dimensions. We 



then apply the appropriate Kogut-Susskind phases [38(, the exact formulation of this phase 
factor depending on the gamma-matrix conventions used in the simulation code. 

Disconnected correlators are by nature noisy; they are directly sensitive to the fluctuations 
in the fermionic and gluonic sea. Consequently it is essential to extract more measurements 
from each configuration than can be gleaned from a single point source inversion. In order to 



achieve this we use the stochastic source method [39|, |40|, |41|, |42J . A set of iV src independent 
noise source vectors {r] a } (a = 1,2,... N src ) is chosen and normalised so that it has the 
orthogonality property 

(vfVj)v = Sofia (20) 
where () v denotes the expectation value over all random noise rj and i labels generically the 
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appropriate vector components such as space, Euclidean time, colour etc. Averaging over 
a sufficiently large number of samples iV src , the set of noise vectors then approximates the 
desired expectation value (1201) and so reliable unbiased estimators for operators sums can 
be constructed. On a given configuration, and on each time slice t we calculate 

■let 

where the sum over i is restricted to the subset of lattice points on timeslice t and, from now 
on, () v denotes the expectation value estimated from the average over N STC noise vectors. 

The Kentucky group has found analytically that among real noise sources, Z(2) noise 
sources should offer the minimal variances for determining propagators 39|. We initially 
tested both complex Z{2) and Gaussian noise sources on a small number of 16 3 x 32 Nf = 2 
lattices and found that Gaussian noise sources produced slightly smaller errors than Z(2) 



noise sources (see 43 1). Based on this preliminary result we employed Gaussian volume 
sources for further simulation. 

To investigate somewhat more thoroughly the apparent discrepancy between this and the 
Kentucky group's result we looked at a larger sample (-/V c f g = 67) of 2 + 1-flavour (3 = 6.76 
lattice configurations, and measured the light quark disconnected correlator D qq (At) using 
both Z(2) and Gaussian noise source vectors. Since the size of the errors of disconnected 
correlators is found to be almost independent of the magnitude of the correlator itself, we 
averaged the error over the timeslice separation At allowing us to quote "errors on the 
errors" . In Fig. [2] the top two curves are the errors on the 75 <8> 1 disconnected correlator 
measured with Gaussian noise (solid line) and with Z(2) noise (dashed line) as a function 
of 1/N mc . It is not apparent that either type of noise presents a large advantage over the 
other in the relevant region of iV src > 8. 

These figures display the interesting property, previously noted in (44)], that the errors 
of disconnected correlators measured with the volume noise source method appear to be 



inversely proportional to N aic , rather than to y/N SIC as one might naively assume. In [45] 
Foley et al. attribute this to a property of the dilution method, however we find it to be 
a generic property of volume sources as well. This can best be understood by considering 
that while the loop operator on a given configuration (121]) 

n sic r 



TV src 



N 

1 ' src , 
a=l 



(22) 
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0.0005 



0.0004 



g 0.0003 
— 



^ 0.0002 
— 

0.0001 



Q-^D qq Z(2) 
□ D Gaussian 

D ss Z ( 2 ) 
A- A D Gaussian 




1/N 



FIG. 2: A comparison of Gaussian and Z(2) noise sources for 67 20 3 x 64 Nj = 2 + 1 lattices for 
f3 = 6.76 am = 0.01, 0.05. The plot shows the standard error on the disconnected correlator D(At), 
averaged over At, versus inverse number of noise sources for the standard (75 ® 1) operator. 



has N src terms in the sum over noise sources, the disconnected correlator on that configura- 
tion has iViL terms: 



V l5m (At) = 15m {t)0 1&m (t+ At) 

iVsrc N s 



-j^ -l * arc - 1 ■ arc 

src a=1 b=1 

so the variance likewise decreases like 1/N t 



vf^smMr/n] 



V^^iM-Wt 

.k,let+At 



(23) 



src" 



C. Variance reduction 



Results from the TrinLat collaboration, using Wilson like fermions, indicate that vari- 
ance is reduced by using "diluted" noise sources when using stochastic sources to measure 
connected correlators [45(. As the name implies, such diluted sources are non-zero only on 
some subset of the lattice. It is then important to use a set of sources for which the non-zero 
subsets together span the entire lattice volume. The subsets may, for example, be times- 
lices, colors, or hypercube corners. We tested several noise dilution schemes in the present 
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0.0004 



G — O g5xl Gaussian noise 
B--H g5xl Z2 noise 
G — o VKVR g5xl Gaussian noise 
B--S VKVRg5xl Z2 noise 



Q 
— 

O 

fcl 



60 

— 

> 



0.0002 



0.0001 




src 



FIG. 3: The dependence of the error on D qq {t) with 1/N src for the 67 lattice configurations described 
in Fig. [21 The comparison is for Z(2) and Gaussian noise for both standard and VKVR-improved 
(75 (g) 1) operators (see text). 



application and found no advantage over using traditional volume-filling sources (see 46|). 

In jsj, Venkataraman and Kilcup introduce a variance reduction technique that is appli- 
cable for the the (75 <8> 1) staggered meson operator. The Venkataraman-Kilcup variance 
reduction (VKVR) trick uses the fact that the staggered Dirac operator M = [Jp +m) and 
its Asqtad-improved variants (and their inverses) connect sites separated by an odd number 
of links. Likewise M^M and its inverse connect sites separated by an even number of gauge 
links. Here we recall that the quark and anti-quark are displaced from each other by four 
links by the (75 ® 1) operator. Now 

ri x M~ l ri y 

is clearly 

7] X (M ] M)~ l M^7] y = r] x (M^M)- 1 (-]p+m)ri y . (24) 

However, since the sites denoted by x and y (quark and antiquark) are 4 links apart, the 
term proportional to Jp has zero expectation value so 



( Vx M-\) v = m( Vx {M^M)-\) v . 



(25) 



13 



4e-05 





E 


1--""" 


C^)D qq Z(2) 

□ D Gaussian 

D ss Z ( 2 ) 
A- -A D Gaussian 
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0.2 



0.4 



0.6 



1/N 



0.8 



"T 
_-t] 




FIG. 4: Dependence of the averaged standard error of D qq {t) and D ss {t) on l/iV src for the 67 lattice 
configurations described in Fig.. [2j VKVR-improved (75 ® 1) operators are used. 

It turns out, however, that the right-hand expression has a significantly decreased variance 
with respect to stochastic noise. There is no additional computational cost since we are 
already computing <fi = M _1 ?7, and we can express the right hand side of (1251) as m(0^0). 

As the VKVR trick should be applicable to disconnected loops for any meson operator 
with the quark and antiquark separated by an even number of lattice spacings, for example 
the (1 ® 1) scalar meson. 

The error on elements of the disconnected correlators is composed of a stochastic part 
- which decreases like 1/N STC as mentioned above — and a gauge part, which remains 
after N* rc 



> 00 and decreases as 1/ wN c { g . The effect of the VKVR trick on disconnected 
correlator errors is dramatic. Fig. [3] shows the error on the disconnected correlator (averaged 
over time separations) as a function of the inverse number of noise sources. With just a 
few noise sources the VKVR trick nearly eliminates all of the stochastic component of the 
error on D qq . At N src = 64 the magnitude of the errors on the VKVR and normal 75 <8> 1 
disconnected correlators become comparable (Fig. [3] and Fig. [5]). 

To calculate the connected and disconnected correlators of the 75 <8> 1 state we wrote 
specialised routines for the Chroma software system j^ . 
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10 20 

t 



FIG. 5: A comparison of disconnected correlators for (3 = 6.76 am = 0.01 obtained with and 
without using the VKVR variance reduced operators for iV src = 64. 

On each configuration, for the light and strange quark masses, we calculated the connected 
correlator for the singlet pion (75 ® 1). We also calculated the mean of Tr [A 75(g) iM _1 ] on 
each timeslice using N STC = 64 Gaussian noise sources for each. The latter we use to calculate 
disconnected correlators. We then use all of these quantities to make the D/C ratio as in 
Equation [16j We fit R(t) to the functional form in ([9]) to determine the splitting between 
the pion and singlet masses. 

IV. RESULTS 

A. Quenched analysis 

With the availability of improved algorithms and significantly more powerful computa- 
tional resources, the quenched approximation to lattice QCD no longer plays such a sig- 
nificant role. However, it remains instructive to study singlet pseudoscalar mesons from 
quenched lattice QCD calculations. In particular the Witten-Veneziano relation relates the 
topological susceptibility, in the large N limit of the SU(N) pure gauge theory, to a contri- 
bution to the mass of the rj meson. A comparison between quenched and full-QCD systems 
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FIG. 6: Dependence of D qq (t) on N STC for 658 = 6.76 am = 0.01 20 3 x 64 lattices. 
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FIG. 7: Disconnected correlators for 658 j3 = 6.76 am = 0.01 20 3 x 64 lattices. 



highlights the effect of dynamical sea quarks. This is particularly instructive in studying 
flavour singlet mesons. Further, since quenched configurations are cheap to generate, one 
can use very large ensembles to illuminate the nature of the statistical fluctuations of singlet 
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quantities and possible long auto-correlation times associated with slow topological modes. 

We calculated connected and disconnected correlators on the quenched gauge configura- 
tions at (3 = 8.00 (neglecting 8 thermalisation configurations) with a valence mass of 0.05. 
In Figure [8] the ratio R(t) is plotted for all the quenched configurations, as well as for subsets 
of the quenched configurations. Quenched staggered QCD involves no determinant fourth- 
roots so there are no obvious theoretical reasons for a significant deviation from (TlOT) for 
the MILC ensemble, when lattice artifacts and finite size effects are neglected. The original 
MILC coarse quenched ensemble {(3 = 8.00) had 408 configurations, and upon extending 
the quenched ensemble as described in Section II III and Table [TTJ, we see from Fig. E] that the 
D/C ratio calculated with the full 6154 configurations is more linear and, as such, in better 
agreement with the expression ffTUl) for At < 15. 
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FIG. 8: D/C ratio for 6154 (3 = 8.00 quenched lattices with valence quark mass am = 0.05, and for 
11 subsets of 400 configurations. We highlight five of these subsets, including the original MILC 
configurations, as disagreeing with the mean by more than one a. 



In Fig. [9] we show the effective mass plots for the light 75 ® 1 and 75 ® 75 non-singlet 
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FIG. 9: Effective masses for the 75 (g> 75 and 75 (g> 1 pseudoscalar mesons for the quenched data 
set (/3 = 8.0). 

pseudoscalar mesons. We computed correlators for the 75 ® 75 pseudoscalar mesons only on 
the original data set, but computed the 75 <g> 1 connected and disconnected correlators on 
the full extended data set. 

In Table II III we report the masses for the 75 ® 1 and 75 <S> 75 pseudoscalar mesons. The 
fits included the full correlation matrix in time and we report both one- and two-cosh mass 
fits. The MILC collaboration obtained the mass 0.46043(20) for the mass of the 75 ® 75 
pseudoscalar meson, using Coulomb gauge fixed wall sources [32] . This is in good agreement 
with our value in Table 1 1 1 1 L As expected, the inclusion of the second state allows us to fit 
much closer to the origin [141 ]. 

Quenched chiral perturbation theory [sl,!^ makes predictions for the ratio of disconnected 



and connected correlators: 



= {ml-am 2 NP ) t + m 2 + am 2 NP 
2m NP 2m 2 NP 



where a is the parameter of the kinetic term of singlet pseudoscalar meson [48| , and mo is 
the difference between the masses msp and m^p. Comparing (1261) with (flUj) we see that it 
defines the parameters A' and B' in fflOl) . 
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TABLE III: Masses for the light pseudoscalars mesons (using connected correlators only) from the 
quenched data set ((3 = 8.00). The 75 (g) 1 correlators were measured with 8 times the statistics of 
the 75 <g> 75 channel. 



Channel 


t-region 


ami 


arri2 


X 2 /dof 


75 ® 75 


13-26 


0.4601(6) 




16.8/ 12 


75 ® 75 


5-25 


0.4603(6) 


1.88(56) 


23.1/17 


75 ® 1 


11-24 


0.5182(3) 




14.2 / 12 


75 ® 1 


5-23 


0.5180(3) 


1.27(5) 


17.1/16 



he value of a is expected to be small, so we neglect it 



indeed Bardeen et 



al. |49|] estimate that a = 0.03 ± 0.03 and advocate setting a = 0. 

Fitting the quenched D(t)/C(t) data to the model in ([TO]) , we obtain a slope of B' 
0.068(3). Using the mass of the 75 ® 1 pseudoscalar meson from Table IIHI and B' — ^ 
we obtain m = 0.76(2) GeV, where we have multiplied by y/nj = \/3 



In Fig. [TUlour value for mo and the values from Bardeen et al. 49j, |51[ are plo t 



2 — • 

NP 



50 



49 



;ed against 
at = 



the square of the pion mass in units of r®. The results from Bardeen et al. 
5.9 were used with the value of ro/a from the ALPHA collaboration 52]. 

Estimates for the value of mo, in the chiral limit, can also be obtained from the Witten 
Veneziano relation 



m n 



4N f XT 
f 2 



where the topological susceptibility xt is evaluated in the 
theory and using /„. with a normalisation of 132 MeV. Shore 



(27) 

arge N c limit of the pure gauge 
53| | describes a modern approach 



;o the use of the topological susceptibility in the study of r] and 7/ mesons. Lucini and Teper 



54i | have shown that the topological susceptibility is only weakly dependent on the number 



of colours. 



The MILC collaboration [55j obtained XtTq = 0.0543(28) at /3 = 8.00 , and xn 



0.0569(26) at (3 — 8.40. These numbers are in good agreement with other quenched results 
such as those by Del Debbio et al. [56], who obtain Xrrt — 0.059(3) in the continuum limit. 
There has also been a recent calculation by Durr et al. 57] . who obtain Xr r o — 0.0524(7) (6) 
in the continuum limit. 
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FIG. 10: Mass dependence of the mo from Bardeen et al. 49f] and this calculation in units of ro- 
There are no factors of n/ included. 

Using XT r o = 0.0543(28) in (1271) with Nf = 3, and taking the value of r = 0.467 fm from 
the HPQCD and MILC collaborations, one obtains m\ = (1.09 GeV) 2 . Using the value of 
r = 0.5 fm, instead, one would obtain = (.952 GeV) 2 . 

The value of m can also be estimated from 19], assuming that the rj is a pure SU(3) 
octet and the rj' is the SU(3) singlet. 



77V + m n — 2m 



K 



(28) 



This method gives a value of wig = (0.854 MeV) 5 



The mass dependence of the mo parameter has recently been discussed by Sharpe 58]. 



In quenched lattice calculations, Kuramashi et al. 



501 ] found that the m parameter had a 



negative slope with respect to the quark mass where it is more usual to see a positive slope 



see the results by MILC 



14J). Sharpe's predictions were tested by Bardeen et al. [5l|. Liu 



and Lagae 59J found that a novel mass-dependent renormalisation factor removed the bulk 

□ 

of the mass dependence of tuq observed in the study by Kuramashi et al. [501 ] . 

Our result mo = 0.76(2) GeV described above, with the valence quark mass set approxi- 
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mately to that of the strange quark, is consistent with other estimates. A chiral extrapolation 
would be required to make a more precise comparison. 



B. Unquenched analysis 

In Fig, mi we present effective mass plots for the (75 (g> 1) correlators, with the connected 
and sum of connected and disconnected correlators, for the (3=6.76, m q /m s = 0.01/0.05 
data set. The effective mass from a correlator of a disconnected loop at the light quark mass 
with a disconnected loop at the strange quark mass is also included in Fig. [TTJ 

In Table IIVI we report single- and double-exponential fits to the connected 75 (g> 1 pseu- 
doscalar correlators. MILC obtained the masses of the 75 ® 75 pseudoscalar mesons to be 



0.22446(22) and 0.49443(25) in lattice units for the light and strange quarks respectively 32]. 
At this lattice spacing the light 75 ® 1 pseudoscalar meson is split from the 75 ® 75 pseu- 
doscalar meson by about 200 MeV. This mass splitting is caused by taste violations and is 
expected to go to zero as the continuum limit is taken. Indeed, the MILC collaboration has 
presented evidence that the mass splittings between different pseudoscalar mesons go like 
0(a 2 a 2 ), as expected, for this action jl^ . 

If the mass splitting between the connected 75 ® 1 pseudoscalar mesons and the 75 ® 75 
correlator was taken as a systematic error at this lattice spacing, then this would imply the 
error in the 77 mass at this lattice spacing is roughly 30%. 

In the twisted mass formalism there is a flavour symmetry breaking term that causes 
a mass splitting between the n° and 7i~ + mesons at non-zero lattice spacing. Although the 



physics of flavour symmetry breaking is different for staggered fermions [36j and the twisted 
formalism, it is encouraging that experience with the twisted mass formalism has shown 
that the disconnected diagrams reduce the mass splitting between the masses of the ir° and 
7r + mesons, caused by isospin violation outside the continuum limit, over the estimate from 



the connected correlators 



60 1 



We see from Table |V] that consistent ground states are determined from the S75S and 
<?75<? operators when the disconnected diagrams are included. This is non-trivial, partic- 
ularly because the masses from the connected correlators are very different for light and 
strange quarks. This is true because the physical rj contains both sj^s and 5755 quark 
content, so either interpolating operator should couple to the r\ meson and produce the 
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ations enforced a consistent mass for the ss and qq operators using 
I2I ] . The contribution from the disconnected part of the correlators 



same ground state mass. The 77 — rj' mixing is further discussed in Section [TTJ The recent 
CP-PACS/JLQCD calcul 
their variational method 
appears to compensate for the increase of the masses from the connected part. 

The correlator between a light and strange pseudoscalar loop should also couple to the 
r) and rj' mesons. Fig. [11] shows that the correlator between a light quark loop and strange 
quark loop (off-diagonal part of (TT81) ) seems to have a plateau at a mass different to the value 
from light or strange correlators (diagonal parts of (|T8l) ). It is not clear why this happening 
- Aubin and Bernard have published propagators for the neutral rj from staggered chiral 
perturbation theory 6l|], but unfortunately these do not include terms for 77 — rf mixing 
that are essential for this analysis. One obvious possibility is that the ground state is 
not being isolated for the singlet correlators in Fig. [TTJ It has been noted before by the 
SESAM collaboration and Bardeen et al. 49|, [5l| that the majority of the excited state 



contamination in a flavour singlet pseudoscalar correlator is from the connected part of the 
correlator. Hence the correlator of a strange loop with a light quark loop could approach the 
ground state at a different rate to correlators with a connected contribution. Given that the 
signal for correlators that include a disconnected contribution in Fig. [TTJ dies in the noise 
beyond timeslice 7, it is difficult to tell where the true plateau lies. 

We tried the variational analysis discussed in Section III Bl The y^jdoj of the fits were 
larger than 7. The problem is, as we discussed above, that the correlators for a light 
loop to strange loop seem to be plateauing at different mass to the singlet light and strange 
loops. The variational fitting formula (IT9l) assumes that each element of the smearing matrix 
couples to the same ground state. We also tried looking at eigenvalues of the matrix, but the 
two eigenvalues just reproduced the diagonal and off diagonal correlators. One possibility 
would be that there is a normalisation problem between the disconnected and connected 
loops. However the magnitude of the light disconnected loop is much larger than that of the 
strange disconnected loop, so a consistent change in the normalisation of the disconnected 
loops makes it hard to get the effective mass of the light-light correlator to plateau at the 
same mass as the strange-strange correlator. This is also a prerequisite for the variational 
analysis. 

The analysis of the non-singlet an correlator using improved staggered fermions turned 
out to be possible, but non-trivial 62J, |63|]. In this case the ground state for the singlet 



22 



0.8 - 



0.6 



0.4 



0.2 



* 



St 



21 



[] 



if 

□ 



[] □ 



« II <> *> 



o 


m = 


0.01 conn only 


□ 


m = 


0.01 conn + disc 




m = 


0.05 conn only 


* 


m = 


0.05 conn + disc 


< 


disc 


(m = 0.05) to disc (m = 0.01) 



# # * 



10 



12 



14 



16 



20 



22 



24 



FIG. 11: Effective masses for the (75 ® 1) channel for the light and strange quarks for the f3 = 
6.76, m=0.01/0.05 data set. 

pseudoscalar channel is the 77 that is stable under the strong interactions. 

The numbers for the lightest mass in Table |V] are about 0.6 in lattice units with weak 
quark mass dependence. This corresponds to a mass in physical units of 990 MeV — the 
ground state in this channel should be the 77 with a mass of 548 MeV. The masses in Table IVl 
are probably contaminated by excited state contributions, so are only an upper limit. A fit 
to the mixed disconnected correlator of a light quark disconnected loop at am q = 0.01 with a 
strange quark disconnected loop at mass am s = 0.05 in Fig.[Tl]gives a mass around 600 MeV. 
Unfortunately because this number does not agree with that from the diagonal correlators 
we do not use it to quote a physical mass for the lightest flavour singlet pseudoscalar meson. 

Fig. [12] shows the ratios of disconnected to connected correlators as defined by ffl6|) in 
the non-degenerate unquenched case (the quenched data is also shown for comparison). The 
ratio plots show that the errors get very large as t gets large. This makes the test that the 
D(t)/C(t) tends to 1 in the large At limit difficult to judge from looking at the graphs. With 
the unquenched data we constructed both R(t)su3 (defined by (|T6l)). and the single-flavour 
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TABLE IV: Masses for the light pseudoscalars mesons for the unquenched data sets obtained from 
(75 ® 1) connected correlators. 



p 


mi/m s 


m v 


region 


ami 


arri2 


X 2 /dof 


6.76 


0.01/0.05 


0.05 


10-17 


0.572(1) 


- 


5.5/6 


6.76 


0.01/0.05 


0.05 


5-11 


0.571(2) 


1.27(11) 


3.5/3 


6.76 


0.007/0.05 


0.05 


11-16 


0.566(2) 


- 


2.6/4 


6.76 


0.007/0.05 


0.05 


5-11 


0.567(3) 


1.19(13) 


4.6/3 


6.76 


0.01/0.05 


0.01 


11-16 


0.363(4) 




4.8/4 


6.76 


0.01/0.05 


0.01 


4-11 


0.358(5) 


1.00(12) 


4.3/4 


6.76 


0.007/0.05 


0.007 


11-16 


0.331(6) 




0.6/4 


6.76 


0.007/0.05 


0.007 


5-11 


0.33(1) 


0.91(34) 


2.9/3 


6.85 


0.05/0.05 


0.05 


12-17 


0.553(2) 




3.8/4 


6.85 


0.05/0.05 


0.05 


4-14 


0.554(2) 


1.33(8) 


8.7/7 



TABLE V: Masses for the light (75 £g> 1) pseudoscalars mesons for the unquenched data sets, from 
the correlators that includes the disconnected diagram. 



p 


m q /m s 


m v 


region 


ami 


X 2 /dof 


6.85 


0.05/0.05 


0.05 


7-11 


0.64(16) 


1.6/3 


6.76 


0.01/0.05 


0.05 


5-11 


0.59(2) 


0.4/3 


6.76 


0.007/0.05 


0.05 


5-9 


0.59(2) 


1.4/4 


6.76 


0.01/0.05 


0.01 


4-9 


0.58(6) 


0.6/3 


6.76 


0.007/0.05 


0.007 


5-9 


0.64(9) 


1.4/3 



R(t) = D(t)/C(t) and fitted both to the model 

R(t) = A + B e' {mA) , (29) 

where m r is the mass difference between the taste singlet pseudoscalar meson and the con- 
nected pseudoscalar correlator. The results for the three fit parameters are given in Table IVT1 
In the table we compare m r , obtained from the ratio, to amf lS . We define the latter to be 
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FIG. 12: D/C ratio for coarse MILC ensembles. For the Nf = 3 and Nf = 2 + 1 ensembles Rsu3 
is displayed. For the quenched configurations the single-flavour R = D/C is plotted. 

the mass of the taste singlet pseudoscalar meson, obtained from the matrix propagator ffTB"]) 
and listed in Table [V] with the mass from the pure connected state (Table HVT) subtracted 
off. 

Within the large 30% errors A is 1 for the light quark as expected from theory. Also the 
values of am^ lS are consistent with am r . Given the noise in the calculation the starting times 
in the fit could only be moved 1 or 2 time units forwards or backwards from those quoted in 
Table IVIj so a more convincing test would require much higher statistics. Unfortunately we 
were unable to get stable fits of the D/C ratio for the Nf — 3 /3 — 6.85 data set, probably 
because the statistics at /3 = 6.85 were lower than for the other /3 values. 



C. Topological Charge 

Smit and Vink 



64| have pointed out that from the Atiyah-Singer index theorem 65j it 
follows that the 75 <S> 1 quark loop is a measure of the topological charge density. Having 
measured this quantity as part of the disconnected correlator calculation we are in a position 
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TABLE VI: Parameters from the fit of (|29|) to the unquenched data. The result for m^ lff is from 
the difference between the mass of the taste singlet and the connected pseudoscalar channel. The 
SU3 entry corresponds to the use of the R ratio in eq. [T5j 



p 


m q /m s 


Mass 


region 


A 


Bo 


am r 




X 2 /dof 


6.76 


0.01/0.05 


0.01 


5 - 13 


1.03(35) 


-0.97(23) 


0.15(11) 


0.22(6) 


6.1/6 


6.76 


0.01/0.05 


SU3 


4 - 11 


1.01(19) 


-1.15(9) 


0.24(9) 


0.22(6) 


6.3/5 


6.76 


0.007/0.05 


0.007 


4- 9 


1.14(39) 


-1.47(23) 


0.24(18) 


0.31(9) 


1.5/3 


6.76 


0.007/0.05 


SU3 


4 - 8 


1.15(30) 


-2.09(91) 


0.35(24) 


0.31(9) 


1.2/2 



to compare the fermionic definition of topological charge 



Q = m Kp (Tr(7 5 M- 1 )) ; 



where renormalisation factor, with the traditional gluonic definition 



(30) 



Q 



64tt 2 



dxe^F^(x)F; a (x) 



(31) 



Alles et al. have performed a similar comparison with Wilson fermions 66] . The MILC 
collaboration has extensively studied the topology on these lattices using gluonic definitions 
of the topological charge |55l. l67|. 
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to measure the topological charge. The 
efore the topological charge 



We used the publicly available MILC code 
gauge fields were cooled by hypercubic blocking [55|, 
was measured. 

A comparison of the two definitions of topological charge on 198 Nf — 2 + 1, (3 — 6.76 
lattices is shown in Fig. [131 an d shows that the gluonic definition of the topological charge is 
strongly correlated with the fermionic definition. This gives us confidence that the staggered 
(Tr^M -1 )) operator, which is a key building block in the SP correlator, is behaving as 



expected. A more empirical comparison between the fermionic and gluonic topologica 



charge 



requires knowledge of the renormalisation factors. Methods have been developed [57J, [7l| to 
determine the matching renormalisation factor between the two definitions of the topological 
charge, based on their assumed equality. 



26 



80 
60 
40 

£ 

§ 

53 



-20 - 

-40 - 

-60 

-80_ 
-30 



o° 



o 
o 



oo 
o 



o 

O O' 



O ' 

o oo° 
o o 

o 



o 
o O 



o 

O u ° CD 

o o °oo 

°°o@ °8 

°<9 o«c? oi° 

O O O OO OO o 



o 



cP ©jo 5 °°° o o 

©j CO 8 o 
Oj o 

o o 



o u o 
o 
>o 



o 



> o 
o 
o 

O O o 



O 

o 

o 
o 

CO 



^o 
o 



o o 



-20 



-10 





Q 



10 



o 

o 

CP 



20 



30 



FIG. 13: Scatterplot of fermionic and gluonic definitions of topological charge on Nf = 2 + 1-flavour 
j3 = 6.76 lattices. Light quark (am = 0.01) loops are used to calculate Tr ^M -1 ). 

V. STATISTICS OF DISCONNECTED CORRELATORS 

Lattice QCD calculations that compute flavour singlet quantities require many more 
gauge configurations than those that only focus on flavour non-singlet quantities. For ex- 
ample, Fig. [8] shows that high statistics are required even in quenched QCD to accurately 
compute the D/C ratio. Although recently there have been impressive algorithmic improve- 
ments in performing unquenched lattice QCD calculations, an increase in the cost of an 
unquenched calculation by a factor of 10 is significant. This motivates a thorough study of 
the statistics of disconnected correlators. 

In Fig. [TU we show a histogram of 393856 measurements (64 timeslices on 6154 quenched 
lattices) of the timeslice loop operator 

L(t) = O ysm (t) = [^ M ~%- ( 32 ) 

igtimcslice t 
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This is for the quenched data set at (3 — 8.00 described in Table HT1 with valence quark 
mass am = 0.05, and shows a clear Gaussian distribution for L. 

However the distribution of the disconnected correlators D(At) constructed from these 
measurements is not Gaussian and standard error estimation requires some care. It is 
straightforward to describe the distribution of D(At) in the limiting case of At = 0. In this 
case (At = 0), D(At) is just given by the square of the loop operator L(t). If independent 
random variables Xi (i.e. L) are Gaussian distributed: 



P(x) 



then Wi = x 2 have a distribution Q{w): 



A exp - 



x 



Q(w) = P(x)^- 
aw 



A 



w 



w 

oxp ( - — 



(33) 



(34) 



Note that the domain of Q(w) is w > 0. This is consistent with the observed distribution 
of quenched disconnected correlator measurements D at At = in Fig. [15j 

In the region of interest, < At < At max (with At max w 15), L(t) and L(t + At) are 
correlated non-trivially and the mean D(At) > 0. In Fig. [161 t ne distribution of D(At = 2) 
is clearly seen to have asymmetric tails. 

Consider two Gaussian-correlated variables x and y with the same variance a 2 . They 
have the joint probability distribution 



P(x, y)=K exp [- (Ax 2 + Ay 2 + 2Bxy)] . 



(35) 



Changing variables to Z{ = Xiyi and 9i = arctan(yj/xj) we get the joint probability for z and 
9: 

dxdy 



Q(z,6) = P(x,y) 



k exp 



dzdB 
2Az 



+ 2Bz 



1 



sin 29 



(36) 



^sin(2^) 

Integrating out 9, treating positive and negative z regions separately (for each it suffices to 
integrate from zero to ±7r/4 and multiply by 4), we get 



Q{z) 



d9Q(z, 



8k exp (-2Bz) 



tx/A 



d9exp 



2Az 
sin(20) 



sin(2^) 



4k exp (-252)^(2^1^1 



(37) 
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where Kq(z) is a modified Bessel function of the second kind. This agrees nicely with, for 
example, D(At = 2) measurements for the quenched f3 = 8.00, am = 0.05 ensemble data, 
shown in Fig. (THl 

We see that in the case of At — > large, D(At) — > and B — > 0. Equation f[3"5"j) then 
factorises into separate x and y parts and 

dxdy 



Q(z,6) = P(x)P(y) 
As expected, we get 



dzdd 



(38) 



Q(z)=4kK (2A\z\). (39) 

This is consistent with the measured quenched disconnected correlators at large time sep- 
aration. Fig. [T7] shows these correlators for At = 20 together with a plot of the modified 
Bessel function form ( 1391) . 

In all cases the most likely measurement is D(At) = 0, irrespective of the mean. Clearly 
the signal we are trying to resolve in taking the mean over disconnected correlator measure- 
ments comes from the asymmetry in the distribution, which is induced by the exponential 
factor in (|37|) . When the mean disconnected correlator is small the distribution is almost 
symmetric and a large proportion of the signal comes from the tails of the distribution. The 
number of data points in the tails of the distribution many standard deviations from the 
mean is far greater than in a Gaussian distribution (Table IVII|) and act as a great lever arm 
on the mean itself. When At is large, both D(At) and C(At) are small and consequently 
the D/C ratio experiences large fluctuations. 

All of this is in marked contrast to the statistics of the connected correlator, the his- 
tograms of which are less peaked and have less pronounced tails (Fig. fT9l) . Another dif- 
ference is obvious — up to differences in the asymmetry, the width of the distribution of 
disconnected correlator measurements remains approximately constant with respect to At, 
whereas the connected correlator distributions are of similar shape with the width roughly 
proportional to the mean. 

The non-Gaussian nature of the distribution of disconnected correlator measurements has 
significant consequences for the interpretation of error estimates on these correlators, and 
upon derived quantities. We use bootstrap sampling methods to estimate the error on the 
mean of quantities of interest. Although the distribution of D is non-Gaussian, the central 
limit theorem ensures that the distribution of the mean D itself is Gaussian. Given large 
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enough samples we should therefore be able to safely use the bootstrap resampling method 
to study the error on the mean using the usual confidence interval methods, but caution is 
called for with smaller samples. 

As an example we consider the 658 configurations of the MILC coarse ensemble with 
(3 = 6.76 and sea-quark masses am = 0.01 and 0.05. We measured the disconnected cor- 
relators and calculated the D/C ratio R(At) as described in Section |IV B[ We dropped 12 
configurations with trajectory index less than 100 and separated the remaining 646 into 
six bins of about 108 configurations each. Fig. [20] shows the D/C ratios for the subsets so 
defined, where the errors were calculated by bootstrap using the usual la (~ 68% confidence 
level) definition. We notice that, at the modest time-separation of At = 10, two subsets (the 
first and last) are nearly a standard deviation from the value obtained with all the data, and 
one subset (the third, highlighted with a filled circle symbol) is nearly two standard devia- 
tions from the value obtained with all 646 configurations. We can immediately trace this to 
the disconnected correlators (in Fig. ETJ the subsets show exactly the same discrepancy in 
the light-quark disconnected correlator D qq ). 

An underestimate of autocorrelation times can cause underestimation of errors, leading to 



an apparent disagreement betwwen subsets of measurements. In [55J the MILC collaboration 
report long autocorrelation times for some 28 3 x 96, a ~ O.lfm ensembles, but not in the 
'coarse' 20 3 x 64 ensembles. Like MILC, we found no measurable autocorrelation of the 
topological charge in the coarse (3 = 6.76, am = 0.01,0.05 ensemble. The related, but 
more relevant quantity, the pseudoscalar singlet disconnected correlator is also suitably 
decorrelated. For D qq (At = 10) we get r ac = 1.09(14). The corresponding timeseries is 
shown in Fig. [22l 

So the autocorrelation time is not the cuplrit here. It is clear from inspection that 
subsets of the timeseries where D qq (At) is particularly deviant from the mean correspond 
to those with a relative abundance (or deficit) of spikes or data points falling in the tails of 
distributions such as in Fig. [TBI 1 

Quenched ensembles exhibit the same features and this is particularly clear from Fig. [83 
We now address the validity of the error bars on the ratio plot for the 400-configuration 



1 The points in the timeseries of D qq (At) measurements are of course the average of N t = 64 correlated 
points from a distribution like Fig. 1181 so the tails and peak are somewhat less severe. 
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subsamples. Fig. [8] shows the ratio R = D/C for the first 400 configurations of each 
stream. The standard error for each set is computed by bootstrap, with a bin size of ten 
configurations. Bootstrap and jackknife give consistent error estimates. From our full 6154- 
configuration quenched ensemble we can make 15 bins of 400 configurations (four of the SO 
bins are not shown in Fig. [S]). Of these, eight have a value of R(10) within its own standard 



error of R(10) = 0.28 (the value for all 6154 configurations). Three give a value between one- 
and two- sigma away, and five have a value between two- and three- sigma away. Assuming 
the means are normally distributed (as discussed above) we would expect that only one of 
the fifteen bins would give a value more than 2a from the "true" value. Although fifteen is 
a small number of bins, this distribution leads us to consider that the size of the error on 
these bins may be slightly underestimated. 

If fl35|) is indeed the form of the distribution of loop operator measurements, then the 
above suggests that fitting the histogram of disconnected correlator measurements to (157)1 
may be a useful method of extracting D(At) = (z). From ( |35|) we easily calculate that 

•2 = A - <* 2 > = jAw (40) 

and 

■y 

(xy) = -— -. (41) 

x 1 2 A 2 - B 2 y 1 

From the Gaussian histogram of loop operator measurements one can extract the overall 

normalization k, and the width of the Gaussian puts a constraint on A and B in the form 

of (|40p . such that there remains effectively one parameter to fit to extract (z) with (j4ip . 

The advantage may be that such a fit would depend less on the fluctuating tails of the 

distribution than the average does. We leave evaluation of this method for future work. 

VI. NUMBER OF CONFIGURATIONS NEEDED 

Flavour singlet lattice spectroscopy is difficult mainly because of the difficulty in getting 
precise determinations of the pseudoscalar disconnected correlators — disconnected correla- 
tors are far more sensitive to fluctuations in the sea than connected correlators. Furthermore 
the pseudoscalar quark loop operator is a measure of the topological charge which is known 



to be plagued with long autocorrelation times 



55(]. 
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Distribution: 


Normal 


P[D(At = 10)] 


K (z) 


P[D(At = 20)] 


la 


0.317311 


0.210280 


0.208994 


0.208147 


2(7 


0.045500 


0.061705 


0.061829 


0.061289 


3(7 


0.002700 


0.019172 


0.019639 


0.019606 


4(7 


0.000063 


0.006297 


0.006460 


0.006751 


5(7 


5.73303e-07 


0.002232 


0.002170 


0.002272 


6(7 


1.97318e-09 


0.000779 


0.000740 


0.000779 


7(7 


2.55962e-12 


0.000322 


0.000255 


0.000297 


So- 


1.22125e-15 


0.000145 


0.000088 


0.000104 


ger 




0.000046 


0.000031 


0.000038 


10a 




0.000025 


0.000011 


0.000025 



TABLE VII: The integrated fraction of data occurring within a radius of multiples of the standard 
deviation a from the mean. Columns two and four are analytic normal (Eq. [33]) and modified 
Bessel (Eq. [39]) distributions, respectively. Columns three and five are the observed distribution 
of D(At = 10) and D(At = 10) for 393856 measurements on quenched lattices. 

We treat this work as a first step toward a more systematic exploration of pseudoscalar 
flavour singlets with staggered quarks. Results such as those depicted in Figures [20] and M 
and the difficulty in fitting indicate that far more configurations are necessary. 

On average, over an ensemble, the error on the 75 ® 1 loop operator is independent of 
timeslice t. A consequence of this is that the product 

75 ®i(*)0 7655 i(* + A£), (42) 

and hence the disconnected correlator D(t), has an error that is roughly constant in t while 
the size of the correlator decreases over many orders of magnitude as t increases. See, for 
example, Fig. [7] 

Since the gauge contribution to the statistical error of the disconnected correlator scales 
as 1/ v/iVtfg, it is straightforward to estimate the number of configurations needed to resolve 
the light-light disconnected correlator (for example) out to some time-separation £ to a 
given precision. We use the (3 = 6.76, am = 0.01,0.05 MILC ensemble, which has 658 
configuration separated by 6 trajectories. We have further grouped the lattices into bins of 
10 lattices (thus eliminating 8). So with 3900 trajectories the gauge error on D qq {t) with 

32 



200 



Af(xj=155xexp(-0.0055x 2 ) 




Tr(y 5 ®l) 

FIG. 14: Histogram of 393856 1^75 ® 1 measurements. Analytic curve is an approximation rather 
than a fit. 

our normalization is roughly 10~ 5 , independent of t (see Fig. HJ) In Fig. [23] we use the 
obtained values of D qq (t) to estimate the number of trajectories needed for various values 
of precision e as a function of time separation. The dotted horizontal line represents the 
current 3900 trajectories, showing that we currently have 20% resolution only out to t = 10 
and 10% resolution out to t — 8. 

It is feasible with modern computational resources to produce ensembles using similar 
parameters with ~ 2 x 10 4 trajectories. The position of the dashed horizontal line represents 
the precision one hopes to derive from these extended statistics. We hope that D qq would 
be resolvable to 20% precision to t = 15 and to 10% precision to t — 11. The use of 
anisotropic lattices may be useful in allowing more time slices where the signal is bigger 
than the noise 0, H] ■ However, the tuning of anisotropic lattice actions is non-trivial for 
unquenched lattice QCD calculations. 
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FIG. 15: Histogram of D(At = 0) measurements, and a curve generated with Equation 



The issue of autocorrelation times has been found not to be a problem at the current 
lattice spacing. In the event of future work on finer lattices — a necessary step for controling 
taste-breaking effects — care should be taken. However, the MILC collaboration report 
in 55| that topological charge autocorrelation times grow with fermion mass, so targeting 
lighter quark masses should help ease the problem. We are concerned with correlations of the 
density of TrA 75lgll M -1 between timeslices. The success of MILC's method of subdividing 
large lattices for dealing with slow topological modes suggests that autocorrelation of the 
topological charge density on timeslices is likely to be less of a problem than that of total 
topological charge. 



VII. CONCLUSIONS 



We have investigated a number of different algorithms to compute the disconnected dia- 
grams required for the correlators for the singlet pseudoscalar mesons. We found that the 



34 



4000 



3000 - 



2; 2000 



1000 




-0.001 



0.001 
D(At = 2) 



0.002 



0.003 



0.004 



FIG. 16: Histogram of 393856 D(At = 2) measurements on 6154 quenched lattices, and a curve 
generated with Equation [371 

algorithm proposed by Venkataraman and Kilcup [9] was the most efficient for our purposes. 

Our results for the ratio of disconnected to connected diagrams (Fig. [T2l) do not re- 
ally show a convincing difference between the quenched (fTUj) and unquenched (Q theories. 
At large time separations the error on the ratio becomes large and it is not clear that it 
asymptotes at 1. This is in contrast to the work by Venkataraman and Kilcup 9] where the 
difference between quenched QCD and unquenched QCD for the ratio was clear. It is clear 
that higher statistics are needed for a definite conclusion on the large-time behaviour of the 
ratio of the disconnected to connected correlator, as well as to extract reliable masses. 

We have started to study the effect of the mixing between light and strange interpolat- 
ing operators. This is an essential part of the physics of the r] and rf meson. The older 
unquenched calculations using the clover and Wilson fermion action 74j . that were done 
with quark masses that were heavier than half the strange quark mass, produced results 
consistent with quenched QCD calculations for the majority of quantities. This suggested 
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FIG. 17: Histogram of D(At = 20) measurements, and a plot of a curve generated from Equation 
[39l the theoretical distribution of the product of two completely uncorrelated Gaussian variables. 

that the strange quark would not play a significant role in the dynamics of 2 + 1-flavour 
calculations. However the strange quark clearly plays an important role in the physics of 
the rj and rf mesons. 



The correlators o: 
the quark operator 



the singlet pseudoscalar meson are closely related to eigenvalues of 
751 ] . This calculation has been done at a lattice spacing of 0.12 fm. 
The calculation of the eigenvalues of the Asqtad improved staggered fermion operator in 
quenched QCD at a lattice spacing of 0.09fm, by Follana at al. [76j, do not show convincing 
clustering of the eigenvalues into quartets. However, the clustering of the eigenvalues for the 
HISQ improved staggered action is convincing. This is some evidence that the computation 
of the spectroscopy of flavour singlet pseudoscalar mesons may require gauge configurations 
with finer lattice spacing than used here. 

Although this lattice calculation has not produced the spectacular agreement with exper- 



36 



6000 
5000 



4000 



2; 3000 



2000 - 



1000 - 





-0.003 



-0.002 



-0.001 



0.001 



0.002 



0.003 



D(At = 10) 



FIG. 18: Histogram of 393856 D(At = 10) measurements on 6154 quenched lattices. 

iment that other parts of the improved staggered program have achieved [lI we have 
not yet seen any "show stoppers" . We are currently generating additional configurations to 
increase the statistics and in order to go to a finer lattice spacing. 
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FIG. 19: Histogram of connected singlet pseudoscalar (75 <g> 1) measurements with am = 0.05 on 
6154 f3 = 8.00 quenched lattices. 

http: / / physics.utah.edu/ ~detar /mile, html ) 
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